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Two numerical methods for nonlinear 
constrained quadratic optimal control 
problems using linear B-spline 
functions 
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Abstract 


This paper presents two numerical methods for solving the nonlinear con- 
strained optimal control problems including quadratic performance index. 
The methods are based upon linear B-spline functions. The properties of 
B-spline functions are presented. Two operational matrices of integration 
are introduced for related procedures. These matrices are then utilized to 
reduce the solution of the nonlinear constrained optimal control to a non- 
linear programming one to which existing well-developed algorithms may be 
applied. Illustrative examples are included to demonstrate the validity and 
applicability of the presented techniques. 
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1 Introduction 


Solving an optimal control problem is not easy. Because of the complexity 
of most applications, optimal control problems are most often solved numer- 
ically. Numerical methods for solving optimal control problems date back 
nearly five decades to the 1950s with the work of Bellman [2-4]. Numeri- 
cal methods for solving optimal control problems are divided into two major 
classes: direct methods and indirect methods. 
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In an indirect method, the calculus of variations [14,24] is used to de- 
termine the first-order optimality conditions of the original optimal control 
problem. The indirect approach leads to a multiple-point boundary-value 
problem that is solved to determine candidate optimal trajectories called ex- 
tremals. Each of the computed extremals is then examined to see if it is 
a local minimum, maximum, or a saddle point. Of the locally optimizing 
solutions, the particular extremal with the lowest cost is chosen. 


One of the widely used methods to solve optimal control problems is 
the direct method. There is a large number of research papers that employ 
this method to solve optimal control problems (see for example [5, 6,10, 15, 
25, 26, 28] and the references therein). This method converts the optimal 
control problem into a mathematical programming problem by using either 
the discretization technique [5,6] or the parameterization technique [10, 25, 
26, 28]. 


The discretization technique converts the optimal control problem into a 
nonlinear programming problem with a large number of unknown parameters 
and a large number of constraints [6]. On the other hand, parameterizing 
the control variables [10, 28] requires the integration of the state equations. 
While the simultaneous parameterization of both the state variables and the 
control variables [28] results in a nonlinear programming problem with a large 
number of parameters and a large number of equality constraints. 


In the last several years, various methods have been proposed to solve 
these problems. Yen and Nagurka [32] proposed a method based on the state 
parameterization, using Fourier series, to solve the linear-quadratic optimal 
control problem (with equal number of state variables and control variables) 
subject to state and control inequality constraints. Also Razzaghi and El- 
nagar [30] proposed a method to solve the unconstrained linear-quadratic 
optimal control problem with equal number of state and control variables. 
Their approach is based on using the shifted Legendre polynomials to pa- 
rameterize the derivative of each of the state variables. In [16] Jaddu and 
Shimemura proposed a method to solve the linear-quadratic and the nonlinear 
optimal control problems by using Chebyshev polynomials to parameterize 
some of the state variables, then the remaining state variables and the control 
variables are determined from the state equations. The approach proposed 
in [28] is based on approximating the state variables and control variables 
with hybrid functions. 


In this paper, we present two computational methods for solving nonlin- 
ear constrained quadratic optimal control problems by using linear B-spline 
functions. The methods are based on approximating the state variables and 
the control variables with a semiorthogonal linear B-spline functions [21]. 
Our methods consists of reducing the optimal control problem to a NLP one 
by first expanding the state rate “(t) and the control u(t) as a linear B-spline 
functions with unknown coefficients. These functions are introduced. For 
the approximation of the integral, the operational matrix of integration Ig is 
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given. Two operational matrices of integration are calculated using (i) dual 
basis functions and (ii) interpolation basis functions. 

The paper is organized as follows: In Section 2 we describe the basic for- 
mulation of the linear B-spline functions required for our subsequent devel- 
opment. Section 3 is devoted to the formulation of optimal control problems. 
Section 4 summarizes the application of these methods to the optimal control 
problems, and in Section 5, we report our numerical findings and demonstrate 
the accuracy of the proposed methods. Sections 6 completes this paper with 
a brief conclusion. 


2 Properties of B-spline functions 


2.1 Linear B-spline functions on [0,1] 


The mth-order cardinal B-spline N,,(t) has the knot sequence {...,—1,0,1,...} 
and consists of polynomials of order m (degree m—1) between the knots. Let 
Ni(t) = xX{0,1)(¢) be the characteristic function of [0,1]. Then for each integer 
m > 2, the mth-order cardinal B-spline is defined, inductively by [8,13] 


Nm(t) = (Nm—1 * N,)(t) = i Nm_—1(t = T)N1(7)dr => i Nm_i(t = T)dr. 


(1) 
It can be shown [7] that N,,,(¢) for m > 2 can be achieved using the following 
formula 


Nm(t) = 


m—1 
recursively, and supp[Nin(t)] = [0, mJ. 
The explicit expressions of N2(t) (linear B-spline function) are [7,8, 13] 


t te (0,1), 
No(t) = 2—tteE (1, 2], (2) 
0 elsewhere. 


Suppose Nj.(t) = No(2/t —k),j,k € Z and Bj, = supp[Nj,x,(t)] = clos{t : 
Nj,n(t) 4 O}. It is easy to see that 


By = ([27%k,279(2+k), j,k eZ. 
To use these functions on [0, 1], 
S; ={k: Bj % A [0,1] F 9}, JEEZ. 


It is easy to see that min{S;} = —1 and max{S;} = 2? —1,j € Z. 
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The support of N;x(t) may be out of interval [0,1], we need that these 
functions intrinsically defined on [0,1] so we put 


djn(t) = Nix(O)xoy(), EZ, ke Sj. (3) 


2.2 The function approximation 


Suppose ®;(t) is a (2/ + 1)-vector as 
®;(t) = [45,-1(€), ;0(¢),---,¢j21-1())", 9 € Z. (4) 


For a fixed j = M, a function f(t) € L?[0, 1] may be represented by the linear 
B-spline functions as 


gM T 
fs S- srour(t) = S?Ou(t), (5) 
k=-1 
where 
S= [s-1, 80,---,82”_4]” (6) 
~~ k+1 
se = f(Sar) | eel ene cane (7) 


Note that the functions ¢j,,,(t) satisfy in the relation 


So we have 
@ys(ti) =e, i= om i= -1,...,2@-1, (8) 


where e; is the ith column of unit matrix of order 2” +1 [21]. 


2.3 Two operational matrices of integration 


Suppose 
t 
wh (t) =f &a(r)ar, (9) 
0 


then the integration of vectors ®j,y in (4) can be expressed as 
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o/, =Iym, (10) 


where Ig is (2“ +1) x (2” + 1) operational matrix of integration for the 
linear B-spline functions on [0,1]. We construct Ig using the following two 
methods: 


Method 1. ; 
t= fof @h (dae (11) 
0 


where ® m is the vector of dual basis of ®;, which can be obtained using the 
linear combinations of ¢;,; [22,23] as 


®yy a P-!@,,, (12) 
where 
ce 
173.3 
1 6 3 6 
Ba) Bip OL (ara ON) ee le (13) 
» yo2 4 
oie 
6 3 


Replacing (12) in (11) we get 


Iz = a Hf, (0% (Oat) P-!=K(P~}), (14) 


where 
E= | wo! (t)@7, (t)dt. (15) 
0 


By using Eqs. (9) and (15) we obtain 


E oa 9-(2M+1) 


Sl- e 
RIF ERIS 
IF 


Method 2. 


In this method, we approximate Oo, using linear B-spline functions and 
then construct Iy. Suppose 


Oh = [Li(t) Lo(t) +++ Lamar (t)]” (16) 
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where using Eq. (4) we have 


t 
Lilt) =f dars-a(r)dr, i=1,...2@+1, M eZ. 
0) 


(Ip) ,; = Lal aM I See Bi, pet OM Sy (17) 


where (Is); denotes the ij-th element of matrix I. Final form of this matrix 


is as follows: 


O11---1 
1 Bees 2 

Ty = 2-444) aera (18) 
1 


3 Problem statement 


The problem we are treating is to find the optimal control u*(t) and the 
corresponding optimal state trajectory x*(t) that minimizes the performance 
index 


J = 5x" (tyxlty) +5 f(x HQAHxlt) +4" REu(d) at, (29) 


subject to 
x(t) = £(x(t), u(t), t), (20) 
V(x(to), to, x(ty),t¢) = 0, (21) 
gi(x(t),u(t),t) <0, 1=1,2,...,w, (22) 


where Z and Q(t) are positive semidefinite matrices, R(t) is a positive def- 
inite matrix, t) and ty are known initial and terminal time respectively, 
x(t) = (x,(t))i_, is the state vector, u(t) = (uj(t))F_y is the control vector 
and f,g; (¢ = 1,2,...,w) are nonlinear functions. This problem is defined on 
the time interval t € [to, t7]. Certain numerical techniques (like B-spline func- 
tions) require a fixed time interval, such as [0,1]. The independent variable 


can be mapped to the general interval 7 € [0,1] via the affine transformation 


r= ; (23) 
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Note that this mapping is still valid with free initial and final times. Using 
Eq. (23), this problem can be redefined as follows: 


Minimize the performance index 


J= 5x1 (I)2x(1) + s(t, - wo) | (x"(r)Q(r)x(r) + u" (r)R(r)u(r)) dr, 


24) 
subject to 
= (ty — to)£(x(7), u(r), 73 to, tf), 25) 
We (0), to,x (1 ), tp) = 0, 26) 
gi(x(T),u(7),7;to, tf) <0, *=1,2,...,w, 7 € [0,1]. 27) 
4 The proposed method 
Let 
@y,1(t) = I, ® ®ys(t), (28) 
®y1q(t) = I, ® y(t), (29) 


where J; and I, are | x | and q x q dimensional identity matrices, @j,(t) is 
(2” + 1)-vector, @ denotes Kronecker product [20] and ®y7,(t) and ©,,4(t) 
are matrices of order (2 + 1) x 1 and q(2™ +1) x q. Assume that each 
of x;(t), 7 = 1,2,...,1, and each of u,(t), 7 = 1,2,...,¢, can be written in 
terms of linear B-spline functions as 


Then using Eqs. (28) and (29) we have 


X(t) ~ OF, (t)X, (30) 
u(t) © O4,,,(t)U, (31) 


where X and U are vectors of orders 1(2” + 1) and q(2™ + 1), respectively, 
given by 


Pp Chee, Corres. call 
US|), Us UO, 


Similarly we have 
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x(0) = ir :(t)Ao, (32) 
where Ag is a vector of order 1(2 +1) given by 
T 
Aj= Lar iia eengr | ‘ 


By integrating Eq. (30) from 0 to t we get 
t 
x(t)—x(0) = | Of, 1(7)Xdr ~ (@ PF, (t))(@1Z)X = OF, (IP X, (33) 
0 


where Iy is an operational matrix of integration given in Eq. (14). From 
Eqs. (32) and (33) we obtain 


x(t) © 7; /(t)(Ao + 17 X). (34) 


4.1 The performance index approximation 


By substituting Eqs. (31)-(34) in Eq. (24) we get 


J = (Ao + 1X)" (yy. (1)ZO, (1))(Ao + 15X) 


slty —toy'Ao+ 5)" ( [banal QW Pr l6)d) (Ao + 1%) 


fs atts —to)UF ( | Barg( ORO, (0d) ul. (35) 


Eq. (35) can be computed more efficiently by writing J as 


J =3(Ao + 1X)" (Z@ O(1)®4,(1) (Ao + I5X) 
dr atts — to)(Ao + 15x)? ( | Q(t) ® au ()O% (Oa) (Ay +I2X) 
4 atts = ty)U* (/ R(t) ® by (0) 0% (Oa) U. (36) 


For problems with time-varying performance index, Q(t) and R(t) are 
functions of time and 


i: Q(t) ® jz (t) 4, (t)dt, i R(t) @ Baz (t)&4,(t)dt 


can be evaluated numerically. For time-invariant problems, Q(t) and R(t) 
are constant matrices and can be removed from the integrals. In this case, 
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Eq. (36) can be rewritten as 


J(X,U) =5(Ao +15 X)* (Z@ ®y(1)®4,(1)) (Ao +13 X) 


+ =(ty — to)(Ao +17 X)7 (Q@ P) (Ap + IX) 


+ =(t7 —t))U? (R@P)U. (37) 


Nl Rew] eR 


4.2 The system constraints approximation 


We approximate the system constraints as follows: 
Using Eqs. (30), (31) and (34) the system constraints (25), (26) and (27) 
became 


Oh, (t)X = (ty = to)f (iz (t) (Ao + 17 X), Oh, (t)U, t; to, tr), 38) 


W (GF, 1(0)(Ao + 17X), to, ®h71(1)(Ao + 17 X), ty) = 0, 39) 
@i(®47(t)(Ao + 13 X), OF, 4(HU, tyto, ty) <0, = i= 1,2,...,w. (40) 


We collocate Eqs. (38) and (40) at Newton-cotes nodes tx, 


k-1 
tk = OM 19202961, (41) 


The optimal control problem has now been reduced to a parameter opti- 
mization problem which can be stated as follows: 
Find X and U so that J(X,U) is minimized (or maximized) subject to 
Eq. (39) and 
ira (te) X = (ty — to)f (Pir (te) (Ao + 15%), Of1,q(te)U, th), (42) 
gi( Phy (t i(t k)(Ao + 1} X), Ohy, g(tk)U, tk; to, tr) < <0, (43) 
i=1,2,...,w, &k=1,2,...,2@41. 


Many well-developed nonlinear programming techniques can be used to solve 
this extremum problem (see, e.g. [1,9,11]). 


5 Illustrative examples 


This section is devoted to numerical examples. All problems were pro- 
grammed in MAPLE, running on a Pentium 4, 2.4-GHz PC with 4 GB of 
RAM. Also we solved the obtained nonlinear programming that is minimize 


26 Y. Edrisi-Tabriz, M. Lakestani and A. Heydari 


(or maximize) J(X, U) subject to Eqs. (39), (42) and (43), using ” NLPsolve” 
command in MAPLE program. To illustrate our techniques, we present five 
numerical examples and make a comparison with some of the results in the 
literatures. 

Example 1. This example is adapted from [18]. Find the control vector 
u(t) which minimizes 


c= sf (aj (t) + u(t) dt, (44) 

subject to 
sta] = [oa] [208] + [2] «© a 
Bo ~ ol , (46) 


and subject to the following inequality control constraint 
ju(t)| < 1. (47) 


In Table 1, the minimum of J using the rationalized Haar functions [29], 
hybrid of block-pulse and Legendre polynomials [25], hybrid of block-pulse 
and Bernoulli polynomials [28] and present two methods are listed. In Figure 
1, the control and state variables with the absolute value of constraint’s errors 
for M = 8, are reported. 


Table 1: Estimated values of J for Example 1 


Method J CPUTime 
Rationalized Haar functions [29] 

K=4 8.07473 0.389 
K=8 8.07065 0.546 
Hybrid of block-pulse and Legendre [25] 

N=4,M, =3 8.07059 1.592 
N=4,M, =4 8.07056 4.304 
Hybrid of block-pulse and Bernoulli [28] 

N=4,M=2 8.07058 0.858 
N=4,M=3 8.07055 1.155 
Present method 1 

M=6 8.07056208507474 1.075 
M=7 8.07056206359357 1.453 
M=8 8.07056204949560 1.722 
Present method 2 

M=6 8.07043910532066 0.665 
M=7 8.07053132323846 1.009 
M=8 8.07055438812380 1.341 
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Figure 1: State and control variables and the constraint errors |&1(t) —2(t)| and |@2(t)+ 
22(t) — u(t)| for Example 1 using Method 1 (left) and using Method 2 (right) with M = 8 
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Example 2. Consider the Breakwell problem [12]. The performance index 
to be minimized is given by 


subject to the state equations 
é(t)=22(t), — &a(t) = u(t), (49) 


with the endpoint conditions 


and the state constraint 
xi(t) < 0.1. (51) 


The exact solution to this problem is given by 


2007 20, ¢€ [0,0.3], 
u*(t)=2 0 t € (0.3, 0.7], (52) 


— 2044+ 49, ¢ € (0.7, 1]. 


This example was studied by using pseudospectral method [12] and ChFD 
scheme [27]. Here we applied the proposed method to solve this problem. 
Absolute errors between approximation and exact value of the performance 
index are reported in Table 2. The approximate solutions of 71(t), x2(t) and 
u(t), obtained by linear B-spline functions using method 2 with M = 9 and 
the exact solutions together error bounds |x7}(t) — x1(t)|, |v3(t) — x2(t)| and 
|u*(t) — u(t)| are plotted in Figure 2. This results show that accuracy of our 
method in comparison with ChFD scheme [27] whose result are plotted in 
Figure 3. 


Table 2: Errors of the estimated and exact values of the performance index, |J — J*|, for 
Example 2 
Method 1 Method 2 

M | [J-J*| CPUTime | |J—J*| CPUTime 

6 | 3.67x 10-2 0.053 5.93 x 10-3 0.163 

7 | 1.86x 10-2 0.181 1.48 x 10-3 0.401 

8 | 9.34x 1073 1.034 3.74. x 107-4 1.377 

9 4.68 x 10-3 7.662 9.38 x 10-5 7.400 


Example 3. This example is adapted from [19] and also studied by using 
rationalized Haar approach [26], hybrid of block-pulse and Legendre polyno- 
mials [25], hybrid of block-pulse and Bernoulli polynomials [28] and interpo- 
lating scaling functions [10]. Find the control vector u(t) which minimizes 
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0.08 


t 


Exact 


° 


Approximate 


Exact °* Approximate 
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bs (4) —x,(1) 


0.000016 


3.x 10°6 
2.x 10°% 


0.000014 


0.000012 


=, 0.000010 
7 
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Figure 2: Exact value, approximation of optimal control, state variables and error bounds 


using method 2 for Example 2 with M = 9 
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0.08 
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* 100 | x7 (=x, (0 | 
+ 10 |x" 5()=x,(D | 
| 20" (t)-ue(t) | 
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0.02 


Figure 3: Exact value and approximation errors of |z7(t) — x1(t)| , |v3(t) — x2(t)| and 
|u*(t) — u(t)| using ChFD scheme [27] for Example 2 with M = 35 


d= , (a7 (t) + x3(t) + 0.005u7(t)) de, (53) 
subject to 

ity | = [oa] [2a] + [1] 64 

al ~ | (55) 


w(t) <r(t), (56) 


where 
r(t) = 8(-0.5)? -0.5, O<t<1. 


The computational result for x2(t) using method 2 for M = 8 together 
with r(¢) are given in Fig. 4. In Table 3, we compare the minimum of J using 
the proposed two methods with other solutions in the literature. 


Example 4. We consider the optimal maneuvers of a rigid asymmetric space- 
craft [17]. This example is studied by using quasilinearization and Chebyshev 
polynomials [15] and hybrid of block-pulse and Bernoulli polynomials [28]. 
The system state equations are 
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Table 3: Results for Example 3 

Method J CPUTime 
Rationalized Haar functions [26] 
K = 64,w = 100 0.170115 1.877 
K = 128,w = 100 0.170103 1.983 
Hybrid of block-pulse and Legendre [25] 
N=4,M, =3 0.17013645 0.951 
N=4,M,=4 0.17013640 1.545 
Hybrid of block-pulse and Bernoulli [28] 
N=4,M=3 0.1700305 0.756 
N=4,M=4 0.1700301 0.921 
Interpolating scaling functions [10] 
n=4,r=5 0.16982646 2.251 
n=5,r=5 0.16982636 3.175 
Present method 1 
M=6 0.169672402102247 1.512 
M=7 0.169782602033829 1.607 
M=8 0.169811048165412 1.985 
Present method 2 
M=6 0.170071967582200 0.599 
M=7 0.169885295276034 1.003 
M=8 0.169837051920398 1.141 


u(t) 


Figure 4: Control and state variables and constraint boundary for Example 3 using 


method 2 with M=8 


0.55 


ba(r) =~ Baa(rya(r) + OO, 

-T ug(T) 
#a(1) = Baa (r)aa(r) + 22, 
2 In —I, U (7) 
£3(T) =— L 1(7)@o(T) + 7, , 


ai(r) — (5 x 10-®r? —5 x 10°-*7 + 0.016) < 0, 


21(100) = 22(100) = x3(100) = 0, 
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where JI, = 86.24, Ig = 85.07, [3 = 113.59. The performance index to be 
minimized, starting from the initial states x,(0) = 0.01, x2(0) = 0.005 and 
x3(0) = 0.001 is 


100 
J= : | (ut (7) + u3(7) + u3(7)) dr. 


We use transformation 7 = 100t,0 < t < 1, in order to use our proposed 
method. In Table 4, the results for J using linear B-spline functions, hy- 
brid of block-pulse and Bernoulli polynomials [28] and quasilinearization and 
Chebyshev polynomials [15] are listed. Optimal control and state variables 
and constraint boundary using method 2, for M = 7, are shown in Figure 5. 


Table 4: Results for Example 4 


Method J CPUTime 
Quasilinearization and Chebyshev polynomials [15] 
N=6 0.00536584 0.07 
N=8 0.00534427 0.10 
N =10 0.00534063 0.12 
Quasilinearization and Chebyshev polynomials [15] 

with using 2 subintervals 
M2 = 10 0.00530902 0.36 
Hybrid of block-pulse and Bernoulli [28] 
N=6,M=3 0.00531097 1.89 
N=6,M=4 0.00530263 2.12 
N=6,M=5 0.00530213 2.74 
Present method 1 
M=5 0.00527460682730895 0.55 
M=6 0.00529464663721832 0.67 
M=7 0.00530275422863559 0.71 
Present method 2 
M=5 0.00530817821841613 0.06 
M=6 0.00530851397972318 0.10 
M=7 0.00530847110421464 0.13 


Example 5. Consider the problem of transferring containers from a ship 
to a cargo truck [31]. The container crane is driven by a hoist motor and a 
trolley drive motor. The aim is to minimize the swing during and at the end 
of the transfer. After appropriate normalization, we summarize the problem 
as follows: 


fe 45 f (02(t) + e2(t)) dt 


subject to 
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Figure 5: Control and state variables and constraint boundary for Example 4 using 
method 2 with M=7 


where 


and 


x1 
£2 
x3 
4 
£5 


x6 


t) = 9xa(t), 

t) = 925(t), 

t) = 9z¢6(t), 

t) = 9(ui(t) + 17.2656x3(t)), 

t) = 9ue(t), 

1) = Mult) + 27.07560r3(t) + 205 (t)n6(t) 


X(t) 


(0) = (0, 22,0,0,—1,0]7, 
a(1) = [10, 14, 0, 2.5, 0, 0]7, 


us (t)| < 2.83374, 
— 0.80865 < u2(t) < 0.71265, 


t € (0, 1], 


with continuous state inequality constraints, 


t € (0, 1], 


? 


In Table 5, we compare the solution obtained using the proposed two 
methods with the method of [9] and [28]. 
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Table 5: Results for Example 5 


Method J CPUTime 
Method of [9] 

m=5 0.5366 x 10~? 2.589 
m=7 0.53614 x 10-2 2.607 
m=9 0.53610895 x 10~2 3.002 
m=11 0.5361102700 x 10~? 3.021 
Hybrid of block-pulse and Bernoulli [28] 

N=2,M=2 0.593000 x 10~? 1.904 
N=2,M=3 0.528915 x 10-2 2.125 
N=2,M=4 0.521421 x 10~? 2.305 
N=2,M=5 0.521411 x 10-2 2.663 
Present method 1 

M=5 0.498574175174882 x 10~? 1.815 
M=6 0.503885802644245 x 107? 1.963 
M=7 0.511514185733751 x 107? 2.025 
Present method 2 

M=5 0.516049181578648 x 10—? 1.579 
M=6 0.515021009757565 x 10~? 1.708 
M=7 0.515021009428266 x 10~? 1.869 


6 Conclusion 


In this paper we presented two numerical methods for solving nonlinear con- 
strained quadratic optimal control problems. Two methods are based upon 
the linear B-spline functions. Also several test problems were used to see 
the applicability and efficiency of the method. The obtained results show 
that when the state variables are unknown at the endpoints, then method 
1 is more accurate than method 2 but in all problems method 2 is faster 
than method 1. In total, our two methods are more accurate than existing 
methods in the literature. 
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